Solving a parabolic problem

Let \(\Omega\) be the unit square \([0,1]^2\). Consider the following problem defined over \(\Omega \times [0,T]\):

\[ \begin{cases} \frac{\partial u}{\partial t}-\Delta u = f &in \ \Omega \times [0,T] \\ \quad u = u_0 &in \ \Omega, \ t=0 \\ \quad u = 0 &on \ \partial \Omega,\ t>0 \end{cases} \]

where \(u_0 = \sin( 2\pi x)\sin(2\pi y)\) is the initial condition, \(f(x,y) = 8\pi^2 \sin( 2\pi x)\sin(2\pi y) e^{-t}\) is the forcing term and \(\partial \Omega\) is the boundary of \(\Omega\) where we have prescribed homogeneous Dirichelet boundary condition for all \(t\) grater than \(0\). The exact solution of the previous problem is \(u_{ex}(x,y) = \sin(2\pi x) \sin(2 \pi y) e^{-t}\) whose trend is shown in the following.

## load domain data
data("unit_square", package="femR")

## time steps
t0 = 0.
t_max = 1
dT = 1e-2
M = (t_max - t0)/dT + 1
times = seq(t0,t_max,length=M)

## Spatio-temporal domain
domain = Mesh(unit_square) %X% times

## create Functional Space
Vh <- FunctionSpace(domain)

## define differential operator in its strong formulation
u <- Function(Vh)

## define differential operator in its strong formulation
Lu <- dt(u) -laplace(u)

## forcing term
f <- function(points,times){
  res <- matrix(0, nrow=nrow(points), ncol=length(times))
  for( t in 1:length(times)){
    res[,t] = (8*pi^2 -1) * sin( 2.* pi * points[,1]) * sin(2.*pi* points[,2]) * exp(-times[t])
  }
  return(res) 
}
## Dirichlet BC
dirichletBC <- function(points, times){
  return(matrix(0, nrow=nrow(points), ncol=length(times)))
}

## initial condition 
initialCondition <- function(points){
 return(sin( 2.* pi * points[,1]) * sin(2.*pi* points[,2]))
}
  
## create pde
pde <- Pde(Lu, f)
#> tracemem[0x55b466553610 -> 0x55b468372d48]: Pde Pde eval eval eval_with_user_handlers withVisible withCallingHandlers handle timing_fn evaluate_call <Anonymous> evaluate in_dir in_input_dir eng_r block_exec call_block process_group.block process_group withCallingHandlers withCallingHandlers handle_error process_file <Anonymous> <Anonymous> <Anonymous> <Anonymous> do.call saveRDS withCallingHandlers doTryCatch tryCatchOne tryCatchList doTryCatch tryCatchOne tryCatchList tryCatch

# set initial condition and dirichlet BC
pde$set_initialCondition(initialCondition)
pde$set_dirichletBC(dirichletBC)

## solve problem
pde$solve()

## plot solution
contour(u)